library(Seurat)
library(tictoc)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(cowplot)

load("Colon_0.6_0722.RDS")

levels(immune.combined)

1. '0.IgD+ mature B'
2. '1.IgM+IgD+ mature B'
3. '2.CD4+CTLA4+ T'
4. '3.CD4+ Teff'
5. '4.Myeloid phagocytes'
6. '5.Ly2Z+ Myeloid'
7. '6.IgA+ Plasma B'
8. '7.Gata3+ ILC2,ILC3'
9. '8.RORga+ ILC2,ILC3'
10. '9.Foxp3+ Treg'
11. '10.CD8+ T, NKT'
12. '11.IgD+ mature B'
13. '12.Neu'
14. '13.proliferare B'
15. '14.B plasmablast'
16. '15.CX2CR1+ Mac'
17. '16.Unknown'
18. '17.Unknown'
19. '18.Unknown'

Principle

In this step, we cluster the cells within each cell type again to separate the doublets as outliers. In Some clusters, like cluster 0 and 1, doublets does not form a unique isolated cluster while others do. We remove those doublet minor populations to enhance the doublet removal performance.

cluster 2

tic()
celltype='2.CD4+CTLA4+ T'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='7'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

cluster 3

tic()
celltype='3.CD4+ Teff'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1165
Number of edges: 55853

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.5346
Number of communities: 8
Elapsed time: 0 seconds
4.43 sec elapsed

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='7'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

cluster 4

tic()
celltype='4.Myeloid phagocytes'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1029
Number of edges: 39231

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.6309
Number of communities: 6
Elapsed time: 0 seconds
4.022 sec elapsed

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='4'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

cluster 5

tic()
celltype='5.Ly2Z+ Myeloid'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 756
Number of edges: 27536

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.6967
Number of communities: 6
Elapsed time: 0 seconds
3.17 sec elapsed

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='5'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

cluster 6

tic()
celltype='6.IgA+ Plasma B'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 755
Number of edges: 31786

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.6785
Number of communities: 7
Elapsed time: 0 seconds
3.235 sec elapsed

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='2'|RNA_snn_res.0.8=='6'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

cluster 14

tic()
celltype='14.B plasmablast'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()

options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
    DimPlot(sub, group.by="DF")+ggtitle(celltype),
    FeaturePlot(sub, features="doublet"),
    DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
    ncol=3
)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 240
Number of edges: 9463

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.5096
Number of communities: 4
Elapsed time: 0 seconds
1.916 sec elapsed

doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='3'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"

total

options(repr.plot.width=10, repr.plot.height=8)
DimPlot(immune.combined, group.by="DF")

# Rename cluster id
Idents(immune.combined)<-"SCT_snn_res.0.6"
new.cluster.ids <- c(
"0.IgD+ mature B", 
"1.IgM+IgD+ mature B", 
"2.CD4+CTLA4+ T", 
"3.CD4+ Teff", 
"4.Myeloid phagocytes", 
"5.Ly2Z+ Myeloid", 
"6.IgA+ Plasma B", 
"7.Gata3+ ILC2,ILC3", 
"8.RORga+ ILC2,ILC3",
"9.Foxp3+ Treg", 
"10.CD8+ T, NKT", 
"11.IgD+ mature B", 
"12.Neu",
"13.Proliferating B", 
"14.B plasmablast",
"15.CX3CR1+ Mac", 
"16.Epithelial cells", 
"17.Endothelial cells", 
"18.Fibroblasts")
names(new.cluster.ids) <- levels(immune.combined)
immune.combined <- RenameIdents(immune.combined, new.cluster.ids)
immune.combined[["celltype0627"]]<-Idents(immune.combined)

options(repr.plot.width=10, repr.plot.height=8)
DimPlot(immune.combined, group.by="celltype0627")

# doublet-free 
seu=subset(immune.combined, DF=="Singlet")
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(seu, group.by="celltype0627")

colnames(immune.combined[[]])

1. 'orig.ident'
2. 'nCount_RNA'
3. 'nFeature_RNA'
4. 'percent.mt'
5. 'stim'
6. 'nCount_SCT'
7. 'nFeature_SCT'
8. 'SCT_snn_res.0.6'
9. 'seurat_clusters'
10. 'celltype0627'
11. 'celltype0627_stim'
12. 'bc'
13. 'DF'
14. 'doublet'

Save clean expression profiles for downstream analysis

save(immune.combined, file="Colon_0.6_0730.RDS", compress=T)